
# Sanghyun Hong, February 23, 2019

#rm(list = ls())
###################################
## Simulation Environment Setup  ##
## Simulation Environment Setup  ##
#################################################################
#pathname<-"C:/Temp/SimulationProject/"
hsig_list=c(0,0.125,0.25,0.5,1.0,2.0,4.0);  			
nj_list=c(40, 80);                                    
al1_list=c(0, 1);                       			    
PubBias=50;                 # Bias = 0-100					           
ind=0;                      # 0=indirect, 1=direct          
Simn=simN;
#################################################################
#
#
#
#
#
#################################################################
# Import Packages and set seed
#################################################################
library(metafor)
library(effsize)
library(sqldf) # This pakcage is for SQL
set.seed(3);
#################################################################
#
#
#
#
#
#################################################################
# Load Packages and Set Path
#################################################################
library(splines)
library(stats)
library(ggplot2)
library(gridExtra)
library(cowplot)
pathnameAK<-paste(pathname,"SimCodes/AK_Codes/",sep="")
setwd(pathnameAK)
source("Modified_AllApplications.R")
#################################################################
#
#
#
#
#
###################################
## Primary Study Data Generation ##
## Primary Study Data Generation ##
#################################################################
dgp <-function(al1, sigh,ind,obs){
  x1=runif(obs, min = 100, max = 200);
  x3=x1+rnorm(obs, mean = 0, sd = 50);
  al3=rnorm(1,mean=0,sd=sigh);
  if (ind==0){
    z = 100 + al1*x1 + al3*x3+rnorm(obs,mean=0,sd=100);
  } else if (ind==1) {
    z = 100 + (al1+al3)*x1 + rnorm(obs,mean=0,sd=100); 
  } 
  return (as.data.frame(cbind(z, x1 ,0, matrix(al1+al3, nrow=obs, ncol=1))))
}
###################################
## Meta Analysis Data Generation ##
## Meta Analysis Data Generation ##
#################################################################
dtcollection <- function(al1, sigh, ssize, Bias,ind, obsList){
  output <-  matrix(0, nrow=ssize, ncol=6);
  colnames(output) <- c("id","y","al_se","mj","Significant","popal");
  
  num_publ=ssize*(Bias/100);
  for(i in 1:ssize) {
    output[i,1]=i;
    obs<-obsList[i];
    if (i<=num_publ){
      while (output[i,5]==0){
        data<-dgp(al1,sigh,ind,obs);
        output[i,6]<-mean(data[,4]);
        out <- lm(data[,1] ~ data[,2]);
        
        output[i,2]<-coefficients(out)[2];		
        output[i,3]<-sqrt(diag(vcov(out)))[2];
        output[i,4]<-0;
        output[i,5]<-as.numeric(summary(out)$coefficients[2,3]>1.96);
      }	
    } else if (i>num_publ){
      data<-dgp(al1,sigh,ind,obs);
      output[i,6]<-mean(data[,4]);
      out <- lm(data[,1] ~ data[,2]);
      
      output[i,2]<-coefficients(out)[2];		
      output[i,3]<-sqrt(diag(vcov(out)))[2];
      output[i,4]<-0;
      output[i,5]<-as.numeric(summary(out)$coefficients[2,3]>1.96);
    }
  }
  return(output)
}
#################################################################
#################################################################
#
#
#
#
#
###################################
## Creating Tables for Output    ##
## Creating Tables for Output    ##
#################################################################
matrixRow=length(hsig_list)*length(nj_list)*length(al1_list)*Simn
FE_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(FE_Output)<-c('TrueEffect','m','sigH','PubBias','EstimatedEffect','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
RE_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(RE_Output)<-c('TrueEffect','m','sigH','PubBias','EstimatedEffect','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
WAAP_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(WAAP_Output)<-c('TrueEffect','m','sigH','PubBias','EstimatedEffect','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
WAAP2_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(WAAP2_Output)<-c('TrueEffect','m','sigH','PubBias','EstimatedEffect','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
PETPEESE_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(PETPEESE_Output)<-c('TrueEffect','m','sigH','PubBias','EstimatedEffect','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
AK_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(AK_Output)<-c('TrueEffect','m','sigH','PubBias','EstimatedEffect','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
AK3_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(AK3_Output)<-c('TrueEffect','m','sigH','PubBias','EstimatedEffect','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
SubTable<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=17))
colnames(SubTable)<-c('TrueEffect','m','sigH','FppCnt','WappCnt','WappTotal','WappSample','WappFrac','Wapp2Cnt','Wapp2Total','Wapp2Sample','Wapp2Frac','AK1P','AK3P1','AK3P2','AK3P3','I2')
#################################################################
#################################################################
#
#
#
#
#
###################################
## Simulation Begins Here        ##
## Simulation Begins Here        ##
#################################################################
start.time <- Sys.time();
sim_i<-0;
cnt<-1;
obsList=rep(c(62,125,250,500,1000),100);

for (k in al1_list){ al1=k; set.seed(3);
for (j in nj_list) { m=j; 
for (l in hsig_list){ sigH=l;
for(i in 1:Simn) {
  sim_i<-sim_i+1;
  data=as.data.frame(dtcollection(al1, sigH, m, PubBias, ind, obsList));
  ###################
  ## FIXED EFFECTS ##
  ## FIXED EFFECTS ##
  ###############################################################
  effectFE<-as.matrix(data$y)
  wFE<-1/data$al_se^2
  RegFE<-summary(lm(effectFE~1, weights=wFE))
  u_FE <- as.numeric(RegFE$coefficients[1])
  seFE <- as.numeric(RegFE$coefficients[2])
  p_valueFE <- as.numeric(RegFE$coefficients[4])
  ciFE <- as.numeric(confint(lm(effectFE~1, weights=wFE), level=0.95))
  covFE <- (ciFE[1]<=al1)*(ciFE[2]>=al1)
  FE_Output[sim_i,]<-c(al1,m,sigH,PubBias,u_FE,seFE,u_FE/seFE,p_valueFE, ciFE[1], ciFE[2], 0, 0, 0, u_FE-al1, covFE, mean(effectFE-al1))
  rm(effectFE, wFE, RegFE, u_FE, seFE, p_valueFE, ciFE, covFE)
  ###############################################################
  ###############################################################
  #
  #
  #
  #
  #
  #
  #
  ###################
  ## Random Effect ##
  ## Random Effect ##
  ###############################################################
  effectRE<-as.matrix(data$y)
  wRE<-1/data$al_se^2
  RegRE<-summary(rma.uni(data$y~1,se=data$al_se, intercept=TRUE, data=data, weighted=TRUE, method="DL", level=95, digits=5));
  tau2RE <- as.numeric(RegRE$tau2)
  I2RE <- as.numeric(RegRE$I2)
  H2RE <- as.numeric(RegRE$H2)
  u_RE <- as.numeric(RegRE$b)
  seRE <- as.numeric(RegRE$se)
  p_valueRE <- as.numeric(RegRE$pval)
  ciRE <- as.numeric(c(RegRE$ci.lb, RegRE$ci.ub))
  covRE<-(ciRE[1]<=al1)*(ciRE[2]>=al1)
  RE_Output[sim_i,]<-c(al1,m,sigH,PubBias,u_RE,seRE,u_RE/seRE,p_valueRE, ciRE[1], ciRE[2], tau2RE, H2RE, I2RE, u_RE-al1, covRE, mean(effectRE)-al1)
  rm(effectRE, wRE, RegRE, H2RE, tau2RE, u_RE, seRE, p_valueRE, ciRE, covRE)
  ###############################################################
  #
  #
  #
  #
  #
  #
  #
  ###################
  ## WAAP          ##
  ## WAAP          ##
  ###############################################################
  wappcnt<-1;
  u_WAAP <- as.numeric(summary(lm(data$y~1, weights=1/(data$al_se^2)))$coefficients[1])
  
  APSdata<-subset(data, (abs(u_WAAP/2.8)>=data$al_se))
  if(nrow(APSdata)<2){APSdata<-data; wappcnt<-0;}
  APSdata<-cbind(APSdata, u_WAAP=u_WAAP)
  
  wapptot <- nrow(data);
  wappsam <- nrow(APSdata);
  wappfrac <- wappsam/wapptot;
  
  effectAPS<-as.matrix(APSdata$y)
  wAPS<-1/APSdata$al_se^2
  RegAPS<-summary(lm(effectAPS~1, weights=wAPS))
  u_APS <- as.numeric(RegAPS$coefficients[1])
  seAPS <- as.numeric(RegAPS$coefficients[2])
  p_valueAPS <- as.numeric(RegAPS$coefficients[4])
  ciAPS <- as.numeric(confint(lm(effectAPS~1, weights=wAPS), level=0.95))
  covAPS <- (ciAPS[1]<=al1)*(ciAPS[2]>=al1)
  if(nrow(data)==nrow(APSdata)){APS_Bias<-u_APS-al1;}else{APS_Bias<-u_APS-al1;}
  WAAP_Output[sim_i,]<-c(al1,m,sigH,PubBias,u_APS,seAPS,u_APS/seAPS,p_valueAPS, ciAPS[1], ciAPS[2], 0, 0, 0, APS_Bias, covAPS, APS_Bias)
  rm(u_WAAP, APSdata, effectAPS, wAPS, RegAPS, u_APS, seAPS, p_valueAPS, ciAPS, covAPS, APS_Bias)
  ###############################################################
  #
  #
  #
  #
  #
  #
  #
  ###################
  ## WAAP2         ##
  ## WAAP2         ##
  ###############################################################
  wapp2cnt<-1;
  u_WAAP <- as.numeric(summary(lm(data$y~1+data$al_se, weights=1/(data$al_se^2)))$coefficients[1])
  
  APSdata<-subset(data, (abs(u_WAAP/2.8)>=data$al_se))
  if(nrow(APSdata)<2){APSdata<-data; wapp2cnt<-0;}
  APSdata<-cbind(APSdata, u_WAAP=u_WAAP)
  
  wapp2tot <- nrow(data);
  wapp2sam <- nrow(APSdata);
  wapp2frac <- wapp2sam/wapp2tot;
  
  effectAPS<-as.matrix(APSdata$y)
  wAPS<-1/APSdata$al_se^2
  RegAPS<-summary(lm(effectAPS~1, weights=wAPS))
  u_APS <- as.numeric(RegAPS$coefficients[1])
  seAPS <- as.numeric(RegAPS$coefficients[2])
  p_valueAPS <- as.numeric(RegAPS$coefficients[4])
  ciAPS <- as.numeric(confint(lm(effectAPS~1, weights=wAPS), level=0.95))
  covAPS <- (ciAPS[1]<=al1)*(ciAPS[2]>=al1)
  if(nrow(data)==nrow(APSdata)){APS_Bias<-u_APS-al1;}else{APS_Bias<-u_APS-al1;}
  WAAP2_Output[sim_i,]<-c(al1,m,sigH,PubBias,u_APS,seAPS,u_APS/seAPS,p_valueAPS, ciAPS[1], ciAPS[2], 0, 0, 0, APS_Bias, covAPS, APS_Bias)
  rm(u_WAAP, APSdata, effectAPS, wAPS, RegAPS, u_APS, seAPS, p_valueAPS, ciAPS, covAPS, APS_Bias)
  ###############################################################
  #
  #
  #
  #
  #
  #
  #
  ###################
  ## PET-PEESE     ##
  ## PET-PEESE     ##
  ###############################################################
  fppcnt <- 0
  effectFPP<-as.matrix(data$y)
  wFPP<-1/data$al_se^2
  se1FPP<-data$al_se
  se2FPP<-data$al_se^2
  RegFPP<-summary(lm(effectFPP~1+se1FPP, weights=wFPP))
  ciFPP <- as.numeric(confint(lm(effectFPP~1 + se1FPP, weights=wFPP), level=0.95)[1,1:2])
  if(as.numeric(RegFPP$coefficients[1,3])>1.96){
    fppcnt <- 1
    RegFPP<-summary(lm(effectFPP~1+se2FPP, weights=wFPP))
    ciFPP <- as.numeric(confint(lm(effectFPP~1 + se2FPP, weights=wFPP), level=0.95)[1,1:2])
  }
  
  u_FPP <- as.numeric(RegFPP$coefficients[1,1])
  seFPP <- as.numeric(RegFPP$coefficients[1,2])
  p_valueFPP <- as.numeric(RegFPP$coefficients[1,4])
  covFPP <- (ciFPP[1]<=al1)*(ciFPP[2]>=al1)
  PETPEESE_Output[sim_i,]<-c(al1,m,sigH,PubBias,u_FPP[1],seFPP[1],u_FPP[1]/seFPP[1],p_valueFPP[1], ciFPP[1], ciFPP[2], 0, 0, 0, u_FPP[1]-al1, covFPP,mean(effectFPP)-al1)
  
  rm(effectFPP, wFPP,se1FPP, se2FPP, RegFPP, ciFPP, u_FPP, seFPP, p_valueFPP, covFPP)
  ###############################################################
  #
  #
  #
  #
  #
  #
  #
  ###################
  ## A-K           ##
  ## A-K           ##
  ###############################################################
  AKdata<-as.data.frame(cbind(data$y, data$al_se, data$id))
  attach(setup(AKdata, TRUE), warn.conflicts = FALSE)
  source("EstimatingSelection.R")
  p_valueAK <- 2*pt(abs(Psihat[1]/se_robust[1]), df=(nrow(my.data)-1),  lower.tail= FALSE);
  ciAK<-as.numeric(Psihat[1]+qt(c(.025, .975), df=(nrow(my.data)-1))*se_robust[1])
  covAK<-(ciAK[1]<=al1)*(ciAK[2]>=al1)
  
  AK_Output[sim_i,]<-c(al1,m,sigH,PubBias,Psihat[1],se_robust[1],Psihat[1]/se_robust[1], p_valueAK, ciAK[1], ciAK[2], 0, 0, 0, Psihat[1]-al1, covAK, mean(data$y)-al1)
  Ak1p1<-Psihat[3]
  rm(Psihat, se_robust, p_valueAK, covAK, ciAK)
  detach("setup(AKdata, TRUE)")
  ###############################################################
  #
  #
  #
  #
  #
  #
  #
  ###################
  ## A-K 3         ##
  ## A-K 3         ##
  ###############################################################
  AKdata<-as.data.frame(cbind(data$y, data$al_se, data$id))
  attach(setup(AKdata, FALSE), warn.conflicts = FALSE)
  source("EstimatingSelection.R")
  p_valueAK <- 2*pt(abs(Psihat[1]/se_robust[1]), df=(nrow(my.data)-1),  lower.tail= FALSE);
  ciAK<-as.numeric(Psihat[1]+qt(c(.025, .975), df=(nrow(my.data)-1))*se_robust[1])
  covAK<-(ciAK[1]<=al1)*(ciAK[2]>=al1)
  
  AK3_Output[sim_i,]<-c(al1,m,sigH,PubBias,Psihat[1],se_robust[1],Psihat[1]/se_robust[1], p_valueAK, ciAK[1], ciAK[2], 0, 0, 0, Psihat[1]-al1, covAK, mean(data$y)-al1)
  Ak3p1<-Psihat[3]; Ak3p2<-Psihat[4]; Ak3p3<-Psihat[5];
  rm(Psihat, se_robust, p_valueAK, covAK, ciAK)
  detach("setup(AKdata, FALSE)")
  ###############################################################
  
  SubTable[sim_i,]<-c(al1,m,sigH, fppcnt, wappcnt, wapptot, wappsam, wappfrac, wapp2cnt, wapp2tot, wapp2sam, wapp2frac, Ak1p1, Ak3p1, Ak3p2, Ak3p3, I2RE)
  rm(data)
  cat("Simulation Environment #9: ", sim_i, "/", matrixRow,"\n");
}
cnt <- cnt+1;
}}}
#################################################################
#################################################################
#
#
#
#
#
###############################################
## Making Graphs and Saving Output           ##
## Making Graphs and Saving Output           ##
###############################################################
FE_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueEffect, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from FE_Output group by TrueEffect, m, sigH'))
RE_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueEffect, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from RE_Output group by TrueEffect, m, sigH'))
WAAP_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueEffect, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from WAAP_Output group by TrueEffect, m, sigH'))
WAAP2_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueEffect, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from WAAP2_Output group by TrueEffect, m, sigH'))
PETPEESE_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueEffect, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from PETPEESE_Output group by TrueEffect, m, sigH'))
AK_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueEffect, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from AK_Output group by TrueEffect, m, sigH'))
AK3_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueEffect, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from AK3_Output group by TrueEffect, m, sigH'))

BIAS_FINAL <- cbind(FE_OUTPUT_FINAL$TrueEffect, FE_OUTPUT_FINAL$m, FE_OUTPUT_FINAL$sigH, FE_OUTPUT_FINAL$PubBias,
                    FE_OUTPUT_FINAL$Bias, RE_OUTPUT_FINAL$Bias,
                    WAAP_OUTPUT_FINAL$Bias, WAAP2_OUTPUT_FINAL$Bias, PETPEESE_OUTPUT_FINAL$Bias, AK_OUTPUT_FINAL$Bias, AK3_OUTPUT_FINAL$Bias)
colnames(BIAS_FINAL)<-c("TrueEffect","m","sigH","PubBias","FE_Bias","RE_Bias","WAAP_Bias","WAAP2_Bias","PETPEESE_Bias","AK_Bias","AK3_Bias")


MSE_FINAL <- cbind(FE_OUTPUT_FINAL$TrueEffect, FE_OUTPUT_FINAL$m, FE_OUTPUT_FINAL$sigH, FE_OUTPUT_FINAL$PubBias,
                   FE_OUTPUT_FINAL$MSE, RE_OUTPUT_FINAL$MSE,
                   WAAP_OUTPUT_FINAL$MSE, WAAP2_OUTPUT_FINAL$MSE, PETPEESE_OUTPUT_FINAL$MSE, AK_OUTPUT_FINAL$MSE, AK3_OUTPUT_FINAL$MSE)
colnames(MSE_FINAL)<-c("TrueEffect","m","sigH","PubBias","FE_MSE","RE_MSE","WAAP_MSE","WAAP2_MSE","PETPEESE_MSE","AK_MSE","AK3_MSE")


COV_FINAL <- cbind(FE_OUTPUT_FINAL$TrueEffect, FE_OUTPUT_FINAL$m, FE_OUTPUT_FINAL$sigH, FE_OUTPUT_FINAL$PubBias,
                   FE_OUTPUT_FINAL$coverage, RE_OUTPUT_FINAL$coverage,
                   WAAP_OUTPUT_FINAL$coverage, WAAP2_OUTPUT_FINAL$coverage, PETPEESE_OUTPUT_FINAL$coverage, AK_OUTPUT_FINAL$coverage, AK3_OUTPUT_FINAL$coverage)
colnames(COV_FINAL)<-c("TrueEffect","m","sigH","PubBias","FE_COV","RE_COV","WAAP_COV","WAAP2_COV","PETPEESE_COV","AK_COV","AK3_COV")


SubTable_Ave<-sqldf('select TrueEffect, m, sigH, avg(FppCnt) as FppCnt, avg(WappCnt) as WappCnt, avg(WappTotal) as WappTotal, avg(WappSample) as WappSample, avg(WappFrac) as WappFrac, avg(Wapp2Cnt) as Wapp2Cnt, avg(Wapp2Total) as Wapp2Total, avg(Wapp2Sample) as Wapp2Sample, avg(Wapp2Frac) as Wapp2Frac, avg(AK1P) as AK1P, avg(AK3P1) as AK3P1, avg(AK3P2) as AK3P2, avg(AK3P3) as AK3P3, avg(I2) as I2 from SubTable group by TrueEffect, m, sigH')

write.csv(BIAS_FINAL, paste(pathname,"/Output/SimulationEnvironment09_BIAS.csv", sep=""))
write.csv(MSE_FINAL,  paste(pathname,"/Output/SimulationEnvironment09_MSE.csv", sep=""))
write.csv(COV_FINAL,  paste(pathname,"/Output/SimulationEnvironment09_COV.csv", sep=""))
write.csv(SubTable_Ave, paste(pathname,"/Output/SimulationEnvironment09_Miscellaneous.csv", sep=""))

end.time <- Sys.time();
time.taken <- end.time - start.time;
cat("Total Time Taken: ", time.taken, "\n\n\n")
###############################################################
###############################################################